#!/bin/bash

#set -o pipefail
set -o errexit
export LC_ALL=C

if [ "x$BASE_DATE" == "x" ]; then
    echo "BASE_DATE variable is not initialized" > /dev/stderr
    exit 1
fi

if [ "x$BASE_DATA_PATH" == "x" ]; then
    echo "BASE_DATA_PATH variable is not initialized" > /dev/stderr
    exit 1
fi

# g++ -O3 6_annotation/intersect_genes.cpp -o 6_annotation/intersect_genes
# g++ -O3 6_annotation/intersect_genes_avg.cpp -o 6_annotation/intersect_genes_avg

# build .so
# g++ -O3 -std=c++0x -c 6_annotation/intersect_genes_simple.cpp -fPIC -lstdc++ -o 6_annotation/intersect_genes_simple.o
# g++ -O3 -std=c++0x -shared -o 6_annotation/intersect_genes_simple.so 6_annotation/intersect_genes_simple.o
# rm 6_annotation/intersect_genes_simple.o

function read_and_prepare_input {
    ANNOTATE_DIR=$1
    GENE_ID=$2
    
    if [ $GENE_ID ]; then
        cat <(echo ">"$GENE_ID) /dev/stdin > $ANNOTATE_DIR/gene.faa
    else
        cat /dev/stdin > $ANNOTATE_DIR/gene.faa
    fi
    
    echo -e $(cat $ANNOTATE_DIR/gene.faa | head -n1 | cut -b2-)"\t"$(cat $ANNOTATE_DIR/gene.faa | head -n2 | tail -n1) \
        > $ANNOTATE_DIR/gene.txt
    
    cat $ANNOTATE_DIR/gene.txt | tools/create_bsddb.py $ANNOTATE_DIR/gen.db
}

function find_local_alignments_blast {
    ANNOTATE_DIR=$1

    cat $ANNOTATE_DIR/gene.faa \
        | 2_raw_matrix_blast/run_blast.sh "${BASE_DATA_PATH}/1_databases/${BASE_DATE}/blast_db/GENOMES_DB" \
        | 2_raw_matrix_blast/convert_output.py \
        | awk -F'\t' '$3<0.05' \
        | sort -t$'\t' -k3,3 -g -T . | sort -s -t$'\t' -k1,1 -T . \
        | gzip \
        > $ANNOTATE_DIR/results_0.05_sorted.txt.gz
}

function find_local_alignments_wu_blast {
    ANNOTATE_DIR=$1
    
    cat $ANNOTATE_DIR/gene.faa \
        | 2_raw_matrix_wu_blast/run_blast.sh "${BASE_DATA_PATH}/1_databases/${BASE_DATE}/wu_blast_db/GENOMES_DB" \
        | 2_raw_matrix_wu_blast/convert_output.py \
        | sort -t$'\t' -k3,3 -g -T . | sort -s -t$'\t' -k1,1 -T . \
        | gzip \
        > $ANNOTATE_DIR/results_sorted.txt.gz

    zcat $ANNOTATE_DIR/results_sorted.txt.gz \
        | awk -F'\t' '$3<0.05' \
        | gzip \
        > $ANNOTATE_DIR/results_0.05_sorted.txt.gz
}

function find_global_alignments_nw_extra {
    ANNOTATE_DIR=$1
    TMP_DIR=$1

    mkdir -p $TMP_DIR/nw_input_files
    mkdir -p $TMP_DIR/nw_output_files

    zcat $ANNOTATE_DIR/results_0.05_sorted.txt.gz \
        | 3_raw_matrix_nw/split_genomes.py $TMP_DIR/nw_input_files

    ls $TMP_DIR/nw_input_files | xargs -n1 -P1 -I {} \
        bash -c "zcat $TMP_DIR/nw_input_files/"{}" \
            | 3_raw_matrix_nw/prepare_data.py  \
                    -G ${BASE_DATA_PATH}/1_databases/${BASE_DATE}/genom.db \
                    -g $ANNOTATE_DIR/gen.db   \
            | 3_raw_matrix_nw/nw_aligner_batch \
                    -i ${BASE_DATA_PATH}/3_raw_matrix_nw/simulation/pupy/interpolation.txt \
            | gzip \
            > $TMP_DIR/nw_output_files/{}"

    find $TMP_DIR/nw_output_files/ -type f | xargs -r zcat | gzip > $ANNOTATE_DIR/nw_results.txt.gz

    rm -r $TMP_DIR/nw_input_files
    rm -r $TMP_DIR/nw_output_files
}

function find_global_alignments_nw_inter {
    ANNOTATE_DIR=$1
    TMP_DIR=$1

    mkdir -p $TMP_DIR/nw_input_files
    mkdir -p $TMP_DIR/nw_output_files

    zcat $ANNOTATE_DIR/results_0.05_sorted.txt.gz \
        | 3_raw_matrix_nw/split_genomes.py $TMP_DIR/nw_input_files

    ls $TMP_DIR/nw_input_files | xargs -n1 -P1 -I {} \
        bash -c "zcat $TMP_DIR/nw_input_files/"{}" \
            | 3_raw_matrix_nw/prepare_data.py  \
                    -G ${BASE_DATA_PATH}/1_databases/${BASE_DATE}/genom.db \
                    -g $ANNOTATE_DIR/gen.db   \
            | 3_raw_matrix_nw/nw_aligner_batch \
            | gzip \
            > $TMP_DIR/nw_output_files/{}"

    find $TMP_DIR/nw_output_files/ -type f | xargs -r zcat | gzip > $ANNOTATE_DIR/nw_results.txt.gz

    rm -r $TMP_DIR/nw_input_files
    rm -r $TMP_DIR/nw_output_files
}

function create_vectors_nw {
    ANNOTATE_DIR=$1
    E_DIR=$2
    Z_VALUE=$3
    NORM_MATRIX_DIR=$4
    EXCL_GENOME_IDS=$5
    
    zcat $ANNOTATE_DIR/nw_results.txt.gz | awk -F'\t' '{if ($3>='$Z_VALUE') print $2"\t"$1}' \
        | 4_normed_matrix/apply_genomes_map.py \
            --genomes-map ${NORM_MATRIX_DIR}/map_int_genom.txt \
            --excl-genomes ${EXCL_GENOME_IDS} \
        | 4_normed_matrix/create_matrix.py \
        > $E_DIR/vectors.txt    
}

function find_sim_vectors_raw_old {
    E_DIR=$1
    P_VALUE=$2
    NORM_MATRIX_DIR=$3
    GENOMES_NUM=$4
    
    6_annotation/intersect_genes \
        -m ${NORM_MATRIX_DIR}/matrix_${P_VALUE}_.txt \
        -v ${E_DIR}/vectors.txt \
        -g $GENOMES_NUM \
    | gzip \
    > $E_DIR/joined_raw.txt.gz
}


function find_sim_vectors_raw {
    E_DIR=$1
    P_VALUE=$2
    NORM_MATRIX_DIR=$3
    GENOMES_NUM=$4
    
    6_annotation/intersect_genes_trans.py \
        -m ${NORM_MATRIX_DIR}/matrix_trans_${P_VALUE}_.db \
        -v ${E_DIR}/vectors.txt \
        -g $GENOMES_NUM \
    | gzip \
    > ${E_DIR}/joined_raw.txt.gz
}


function find_sim_vectors {
    E_DIR=$1
    P_VALUE=$2
    NORM_MATRIX_DIR=$3
    
    join -t$'\t' -1 1 -2 2 \
        <(cat ${NORM_MATRIX_DIR}/map_int_gene.txt | sort -t$'\t' -k1,1 -S100M) \
        <(zcat ${E_DIR}/joined_raw.txt.gz | sort -t$'\t' -k2,2 -S100M) \
    | cut -f2,4 \
    | sort -t$'\t' -k1,1 -S100M \
    | gzip \
    > ${E_DIR}/joined.txt.gz

    zcat ${E_DIR}/joined.txt.gz \
        | sort -t$'\t' -k2,2 -S100M -g -r \
        | tail -n1000 \
        | sort -t$'\t' -k1,1 -S100M \
        | gzip \
        > ${E_DIR}/.results_best_1K.txt.gz

    mv ${E_DIR}/.results_best_1K.txt.gz ${E_DIR}/results_best_1K.txt.gz
}

function create_kegg_limited {
    E_DIR=$1
    P_VALUE=$2
    PROB_LIMIT=$3
    
    join -t$'\t' \
        <(zcat ${E_DIR}/joined.txt.gz) \
        ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/gene_kegg_map.txt \
    | awk -F'\t' '{if ($2<'"$PROB_LIMIT"') print $0}' \
    | gzip \
    > ${E_DIR}/joined_kegg_limited.txt.gz

    #rm ${E_DIR}/joined.txt.gz
}

function create_gene_ontology {
    E_DIR=$1
    P_VALUE=$2
    PROB_LIMIT=$3

    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        <(zcat ${E_DIR}/joined.txt.gz | awk -F'\t' '{if ($2<'"$PROB_LIMIT"') print $0}') \
        <(cat ${BASE_DATA_PATH}/1_databases/${BASE_DATE}/genes_go.txt ) \
    | cut -f2,3 \
    | sort -t$'\t' -g -k1,1 -S100M \
    | sort -t$'\t' -s -k2,2 -S200M \
    | awk -F'\t' '{if ($2!=S) {if (S) {print P"\t"C"\t"S}; S=$2; P=$1; C=-1} else {if ($1<P) {P=$1}; C=C-1}}END{if (S) {print P"\t"C"\t"S}}' \
    | sort -t$'\t' -g -k1,1 -k2,2 -S200M \
    | awk 'BEGIN{OFS="\t"; FS="\t"}{$2=-$2; print}' \
    | cut -f1,3 \
    | gzip \
    > ${E_DIR}/.results_gene_ontology.txt.gz

    mv ${E_DIR}/.results_gene_ontology.txt.gz ${E_DIR}/results_gene_ontology.txt.gz
}

function clear_tmp_files {
    E_DIR=$1

    rm -f ${E_DIR}/joined_raw.txt.gz
    rm -f ${E_DIR}/joined.txt.gz
}

function create_gene_ontology2 {
    E_DIR=$1
    P_VALUE=$2
    PROB_LIMIT=$3

    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    zcat ${E_DIR}/results_gene_ontology.txt.gz \
    | sort -t$'\t' -s -k2,2 -S200M \
    | awk -F'\t' '{if ($2!=S) {if (S) {print P"\t"C"\t"S}; S=$2; P=$1; C=-1} else {if ($1<P) {P=$1}; C=C-1}}END{if (S) {print P"\t"C"\t"S}}' \
    | sort -t$'\t' -g -k1,1 -k2,2 -S200M \
    | awk 'BEGIN{OFS="\t"; FS="\t"}{$2=-$2; print}' \
    | cut -f1,3 \
    | gzip \
    > ${E_DIR}/results_gene_ontology2.txt.gz

    mv ${E_DIR}/results_gene_ontology2.txt.gz ${E_DIR}/results_gene_ontology.txt.gz
}

function create_gene_ontology3 {
    E_DIR=$1
    P_VALUE=$2
    PROB_LIMIT=$3

    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        <(zcat ${E_DIR}/joined.txt.gz | awk -F'\t' '{if ($2<'"$PROB_LIMIT"') print $0}') \
        <(cat ${BASE_DATA_PATH}/1_databases/${BASE_DATE}/genes_go.txt ) \
    | cut -f2,3 \
    | sort -t$'\t' -g -k1,1 -S100M \
    | sort -t$'\t' -s -k2,2 -S200M \
    | awk -F'\t' '{if ($2!=S) {if (S) {print P"\t"C"\t"S}; S=$2; P=$1; C=-1} else {if ($1<P) {P=$1}; C=C-1}}END{if (S) {print P"\t"C"\t"S}}' \
    | sort -t$'\t' -g -k1,1 -k2,2 -S200M \
    | awk 'BEGIN{OFS="\t"; FS="\t"}{$2=-$2; print}' \
    | cut -f1,3 \
    | gzip \
    > ${E_DIR}/results_gene_ontology_1e6.txt.gz
}

function create_gene_functions {
    E_DIR=$1
    P_VALUE=$2

    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        <(zcat ${E_DIR}/joined_kegg_limited.txt.gz) \
        <(cat ${BASE_DATA_PATH}/1_databases/${BASE_DATE}/genes_info.txt ) \
    | cut -f2,3,6 \
    | sort -t$'\t' -g -k1,1 -S100M \
    | gzip \
    > ${E_DIR}/results_gene_functions.txt.gz

}


function create_gene_functions_kegg {
    E_DIR=$1
    P_VALUE=$2
   
    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        <(zcat ${E_DIR}/joined_kegg_limited.txt.gz) \
        <(cat ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/gene_functions.txt)  \
       | cut -f2,3,4 \
    | sort -t$'\t' -g -k1,1 -S100M \
    | gzip \
    > ${E_DIR}/results_gene_functions_kegg.txt.gz

}

function create_gene_functions_uniq {
    E_DIR=$1
    P_VALUE=$2

    ##
    ## functions_cnt<TAB>function
    ##
    zcat ${E_DIR}/results_gene_functions.txt.gz \
        | cut -f3 | sort -t$'\t' | uniq -c | sort -k1,1 -nr \
        > ${E_DIR}/results_gene_functions_uniq.txt
}


function create_gene_functions_best_uniq {
    E_DIR=$1
    P_VALUE=$2

    ##
    ## functions_cnt<TAB>function
    ##
    zcat ${E_DIR}/results_gene_functions.txt.gz \
        | head -n1000 | cut -f1,3 | sort | uniq -c | sort -k2,2 -gs | sort -t$'\t' -k2,2 -su | sort -k1,1 -nr | sort -k2,2 -gs \
        > ${E_DIR}/results_gene_functions_best_uniq.txt
}

function create_gene_functions_exp {
    E_DIR=$1
    P_VALUE=$2

    zcat ${E_DIR}/results_gene_functions.txt.gz \
        | cut -f1,3 | sort -t$'\t' -gk1,1 -S200M | sort -t$'\t' -sk2,2 -S200M \
        | awk -F'\t' '{if ($2!=S) {if (S) {print P"\t"C"\t"S}; S=$2; P=$1; C=-1} else {if ($1<P) {P=$1}; C=C-1}}END{if (S) {print P"\t"C"\t"S}}' \
        | sort -t$'\t' -g -k1,1 -k2,2 -S200M | awk 'BEGIN{OFS="\t"; FS="\t"}{$2=-$2; print}' | cut -f1,3 \
        | gzip \
        > ${E_DIR}/results_gene_functions_exp.txt.gz
}

function create_gene_functions_best_uniq_kegg {
    E_DIR=$1
    P_VALUE=$2

    ##
    ## functions_cnt<TAB>function
    ##
    zcat ${E_DIR}/results_gene_functions_kegg.txt.gz \
        | head -n1000 | cut -f1,3 | sort | uniq -c | sort -k2,2 -gs | sort -t$'\t' -k2,2 -su | sort -k1,1 -nr | sort -k2,2 -gs \
        > ${E_DIR}/results_gene_functions_kegg_best_uniq.txt
}

function create_gene_functions_kegg_exp {
    E_DIR=$1
    P_VALUE=$2

    zcat ${E_DIR}/results_gene_functions_kegg.txt.gz \
        | cut -f1,3 | sort -t$'\t' -gk1,1 -S200M | sort -t$'\t' -sk2,2 -S200M \
        | awk -F'\t' '{if ($2!=S) {if (S) {print P"\t"C"\t"S}; S=$2; P=$1; C=-1} else {if ($1<P) {P=$1}; C=C-1}}END{if (S) {print P"\t"C"\t"S}}' \
        | sort -t$'\t' -g -k1,1 -k2,2 -S200M | awk 'BEGIN{OFS="\t"; FS="\t"}{$2=-$2; print}' | cut -f1,3 \
        | gzip \
        > ${E_DIR}/results_gene_functions_kegg_exp.txt.gz
}


function create_gene_pathways {
    E_DIR=$1
    P_VALUE=$2

    ##
    ## probability<TAB><gene_kegg_id<TAB>path_id<TAB>pathname
    ##
    join -t$'\t' -1 2 \
        <(join -t$'\t' -2 2 \
            ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/kegg_pathway_map_full.txt \
            <(zcat ${E_DIR}/results_gene_functions.txt.gz | cut -f1,2 | sort -t$'\t' -k2,2 -S100M) \
        | sort -t$'\t' -k2,2 -S100M \
        ) \
        <(cat ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/pathways_full.txt | cut -f1,4) \
        | awk -F'\t' '{print $3"\t"$1"\t"$2"\t"$4}' \
        | sort -t$'\t' -g -k1,1 -S100M \
        | gzip \
        > ${E_DIR}/results_pathways.txt.gz
}


function create_gene_pathways_avg_nw {
    E_DIR=$1
    Z_VALUE=$2
    GENOMES_NUM=$3

    ## 
    ## probability<TAB>path_id<TAB>intersection<TAB>pathname
    ##    
    join -t$'\t'  \
        <(6_annotation/intersect_genes_avg \
            -m ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/pathway_vectors_avg_wunw_${Z_VALUE}_.txt \
            -v $E_DIR/vectors.txt \
            -g $GENOMES_NUM \
            | cut -f2- | sort -t$'\t' -k1,1 -S100M) \
        <(cat ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/pathways_full.txt | cut -f1,4) \
        | awk -F'\t' '{print $2"\t"$1"\t"$3"\t"$4}' \
        | sort -t$'\t' -g -k1,1 -S100M \
        > $E_DIR/results_pathways_avg.txt
}

function create_gene_pathways_all {
    E_DIR=$1
    PROB_LIMIT=$2

    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' -1 3 -2 1 \
        <(join -t$'\t' \
            <(zcat ${E_DIR}/results_best_1K.txt.gz | awk -F'\t' '{if ($2<'"$PROB_LIMIT"') print $0}') \
            <(cat ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/gene_pathway_map.txt ) \
            | cut -f2,3 \
            | sort -t$'\t' -g -k1,1 -S100M \
            | sort -t$'\t' -s -k2,2 -S200M \
            | awk -F'\t' '{if ($2!=S) {if (S) {print P"\t"C"\t"S}; S=$2; P=$1; C=-1} else {if ($1<P) {P=$1}; C=C-1}}END{if (S) {print P"\t"C"\t"S}}' \
        ) \
        ${BASE_DATA_PATH}/5_pathways/${BASE_DATE}/pathways_full.txt \
    | awk 'BEGIN{OFS="\t"; FS="\t"}{print $2"\t"(-$3)"\t"$1"\t"$4"\t"$6}' \
    | sort -t$'\t' -k2,2nr -k4,4n -S200M \
    | gzip \
    > ${E_DIR}/.results_pathways_all.txt.gz

    mv ${E_DIR}/.results_pathways_all.txt.gz ${E_DIR}/results_pathways_all.txt.gz
}

